
% bananas
% Christopher J. Whalen
% University of Illinois
% 10/11/04 - Date Created
% 10/18/04 - Last Modified

% Purpose - to create a banana brik of the digitized points and including
% the tai fiducials
%           makes bananas green to red depending upon overlap

 function [dig_brik] = bananas(choice,file,xyz,src,det,tai,opt_label,mri_dim,head_mask,brain_mask)     

% fprintf('\nCreating Banana BRIK ...\n');
choice.banana_type = input('\nBanana Type (Comp. Time):\n1 - Discrete (1x)\n2 - Fuzzy (20x)\n3 - Super Fuzzy (34x)\n');
% Calculating ref_banana_factor which is time for computation based on file
% size
if choice.banana_type == 1
    ref_banana_factor = 1;
elseif choice.banana_type == 2
    ref_banana_factor = 20;
elseif choice.banana_type == 3
    ref_banana_factor = 175;
end

%load('C:\Analysis\Matlab\coreg\banana_input.mat')
%save C:\Analysis\Matlab\coreg\banana_input.mat 'tai_origin' 'src' 'det' 'tai_fid' 'opt_label' 'dig_brik' 'xyz'

% define acpc.BRIK size and orientation
dig_brik=zeros(mri_dim.x,mri_dim.y,mri_dim.z);  % RAI

%xyz = round(xyz);
% xyz = xyz - 1;           % Subtracting 1 mm to go back into AFNI coords (?)
% tai.fid = tai.fid - 1;   % Subtracting 1 mm to go back into AFNI coords (?)
% tai.origin=[0 0 0]; % I have no idea why Chris adds this ~line 177

min_dist = input('\nMinimum S-D Distance for Banana (Default Value 25 mm): ');
max_dist = input(  'Maximum S-D Distance for Banana (Default Value 50 mm): ');
if isempty(min_dist)
    min_dist = 25; % minimum allowable source detector distance (e.g. 15 = 15 mm)
end
if isempty(max_dist)
    max_dist = 50; % maximum allowable source detector distance: this is can not be greater than 100 mm
end
mu_a = 0.05; % absorbtion coefficient

% Must eventually change to put in as real parameters as possible
%
n_cube_ref = 100;
n_x_ref = n_cube_ref; n_y_ref = n_cube_ref; n_z_ref = n_cube_ref;
ref_banana = zeros(n_x_ref,n_y_ref,n_z_ref);
s_mrkr_pt = [n_x_ref n_y_ref/2 n_z_ref]; % making a marker point to follow appropriate rotations
banana_list = zeros(n_cube_ref^3,4); % banana_list is x y z coord and intensity value
banana_list(:,4) = 1; % setting the default intensity to be 1
banana_brik = zeros(n_x_ref,n_y_ref,n_z_ref);
banana_intensity_increment = 1; % The increment in the color scheme 4 - Green; 127 - Red
max_intensity = 126;
ref_sd_dist = 100;
banana_radius = ref_sd_dist/2;
max_banana_girth = banana_radius/4;
ref_sd_mid_pt = [n_cube_ref/2 n_cube_ref/2 n_cube_ref];
%banana_lower_bound = banana_radius - max_banana_girth;
%banana_upper_bound = banana_radius + max_banana_girth;
s_ref = [0 n_y_ref/2 n_z_ref];
d_ref = [n_x_ref n_y_ref/2 n_z_ref];
s_ref = s_ref - ref_sd_mid_pt;
s_ref = s_ref/sqrt(s_ref(1)^2+s_ref(2)^2+s_ref(3)^2);
ban_ref = [0 0 -1];

if choice.banana_type == 1 & fopen(file.ref_banana)~=-1
%     fprintf('\nLoading Reference Banana...\n');
    load(file.ref_banana);
elseif choice.banana_type == 2 & fopen(file.ref_banana_fuzzy)~=-1
%     fprintf('\nLoading Reference Fuzzy Banana...\n');
    load(file.ref_banana_fuzzy);
elseif choice.banana_type == 3 & fopen(file.ref_banana_fuzzy_hd)~=-1
%     fprintf('\nLoading Reference Fuzzy Banana...\n');
    load(file.ref_banana_fuzzy_hd);
else
    % Generating Reference Banana
    [banana_list] = gen_ref_banana(choice,n_x_ref,n_y_ref,n_z_ref,banana_radius,max_banana_girth,banana_list,banana_brik);
    if choice.banana_type == 1 
        fprintf('\n Saving Reference Banana ...\n');
        save(file.ref_banana,'banana_list','-mat'); 
    elseif choice.banana_type == 2 
        fprintf('\n Saving Reference Fuzzy Banana ...\n');
        save(file.ref_banana_fuzzy,'banana_list','-mat'); 
    elseif choice.banana_type == 3 
        fprintf('\n Saving Reference High Definition Fuzzy Banana ...\n');
        save(file.ref_banana_fuzzy_hd,'banana_list','-mat'); 
    end
end

% looping through all sources and detectors pairs
%
% fprintf('\nCreating S-D Bananas...\n');
tot_channels = length(src)*length(det);
% fprintf('Total Possible Channels:  %d\n', tot_channels);
% Determining Valid SD pairs to make time estimate
sd_pair = zeros(length(src)*length(det),1);
m = 1; % counter
for a=1:length(src)
    for b=1:length(det)
        for c=1:length(opt_label)
            if strcmp(src{a},opt_label{c})
                s = xyz(c,:);
                %opt_label(k) % line to debug if errors
            end
            if strcmp(det{b},opt_label{c})
                d = xyz(c,:);
                %opt_label(k)  % line to debug if errors
            end
        end
        sd_pair(m) = (d(1)-s(1))^2 + (d(2)-s(2))^2 + (d(3)-s(3))^2;
        m = m + 1;
    end
end
sd_pair = sort(sd_pair);
min_index = find(sd_pair>min_dist^2,1,'first') - 1;
if (max_dist^2)>sd_pair(size(sd_pair,1))
    max_index = size(sd_pair,1);
else
    max_index = find(sd_pair>max_dist^2,1,'first');
end
valid_channels = max_index - min_index; % range of valid channels
percent_valid_channels = 100*valid_channels/tot_channels;
t_comp = valid_channels*0.025*ref_banana_factor;

fprintf('Number accepted channels / total channels:  %d / %d (%-2.1f %%)\n',...
    valid_channels,tot_channels, percent_valid_channels);   
% fprintf('Number Invalid Channels:  %d\n', (tot_channels - valid_channels) );
% fprintf('\nEstimate Time: %-2.1f min\n', t_comp); % Estimating Time
%error_det = cell(tot_channels,7); %trying to detect where errors are in
%rotations

tic
q = 1; % counter for sd valid pairs
for i=1:length(src)  % 8:8 % % This cycles through all possible sources
    for j=1:length(det)  %37:37 % % This cycles through all possible detectors
        for k = 1:length(opt_label)
            if strcmp(src{i},opt_label{k})
                s = xyz(k,:);
                %opt_label(k) % line to debug if errors
            end
            if strcmp(det{j},opt_label{k})
                d = xyz(k,:);
                %opt_label(k)  % line to debug if errors
            end
        end
        sd_dist = sqrt( (d(1)-s(1))^2 + (d(2)-s(2))^2 + (d(3)-s(3))^2);
        sd_mid_pt = [(s(1)+d(1)) (s(2)+d(2)) (s(3)+d(3))]/2;
        if (sd_dist <= max_dist) & (sd_dist >= min_dist)
%             c = rand; % decimating bananas for viewing
%             if c>0.025  %decimating bananas
%                 % do nothing
%             else
            
%             fprintf('\n%d/%d Generating Banana for Source-Detector Pair %s-%s',q, valid_channels, src{i}, det{j});
            banana_brik = zeros(mri_dim.x,mri_dim.y,mri_dim.z); % brik for each particular banana
            q = q + 1;
            % Scaling Down Reference Banana and Removing Extraneous Points
            banana_scaling_factor = sd_dist/n_cube_ref; % scales ref banana to appropriate size
            sd_banana = banana_list*banana_scaling_factor; %defining sd specific banana
            sd_banana(:,4)=sd_banana(:,4)/banana_scaling_factor; % compensates for multiplication in prev step
            sd_banana_mid_pt = ref_sd_mid_pt*banana_scaling_factor;
            sd_banana(:,1:3) = sd_banana(:,1:3) - repmat(sd_banana_mid_pt,length(sd_banana), 1); % making origin equal to midpoint of source detector
            s_mrkr_pt = s_mrkr_pt*banana_scaling_factor; % scaling reference marker
            s_mrkr_pt = s_mrkr_pt - sd_banana_mid_pt;
            p = [0 0 -1]; % Point P gives the vector of the ref source fiber vector
            % Can maybe put this after transform
            % Remove Identical Points in Banana List
            [sd_banana] = rm_identical(sd_banana);
            
            % RBTfm Banana
            % defining three points in ref banana frame
            ref_banana_pts = [0 0 0; s_ref; ban_ref];
            s_vect = -1*(sd_mid_pt-tai.origin)/(sqrt((sd_mid_pt(1)-tai.origin(1))^2+(sd_mid_pt(2)-tai.origin(2))^2+(sd_mid_pt(3)-tai.origin(3))^2));    % direction source fiber is pointing. Now Unit Vector toward tai Origin. this should change if eventually use normal vector
            a= cross(cross(s_vect,(s-d)/sd_dist),(s-d)/sd_dist);
            a = a/(sqrt(a(1)^2+a(2)^2+a(3)^2));
            s_vect_correct = dot(s_vect,a)*a;
            s_vect_correct = s_vect_correct/sqrt(s_vect_correct(1)^2+s_vect_correct(2)^2+s_vect_correct(3)^2);
            dist1 = sqrt((s_vect(1) - s_vect_correct(1))^2 + (s_vect(2) - s_vect_correct(2))^2 + (s_vect(3) - s_vect_correct(3))^2 );
            dist2 = sqrt((s_vect(1) + s_vect_correct(1))^2 + (s_vect(2) + s_vect_correct(2))^2 + (s_vect(3) + s_vect_correct(3))^2 ) ;
            if dist2<dist1 % compensating for going in neg direction
                s_vect_correct = -s_vect_correct;
            end
            % three points in optode frame
            sd_loc = [0 0 0; (s-d)/sd_dist; s_vect_correct];
            

            [R,T] = get_lsrbtfm(ref_banana_pts,sd_loc);
            C = calc_centroid(ref_banana_pts);
            sd_banana(:,1:3) = sd_banana(:,1:3) - repmat(C,length(sd_banana), 1); % placing back in BRIK coord
            %sd_banana_fuzzy(:,1:3) = sd_banana_fuzzy(:,1:3) - repmat(C,length(sd_banana_fuzzy),1);
            for n=1:size(sd_banana,1)
                sd_banana(n,1:3)=(R*sd_banana(n,1:3)' + T)';
            end
            %for n=1:size(sd_banana_fuzzy,1)
            %    sd_banana_fuzzy(n,1:3)=(R*sd_banana_fuzzy(n,1:3)' + T)';
            %end
            %sd_banana_fuzzy(:,1:3) = sd_banana_fuzzy(:,1:3) + repmat(sd_mid_pt+C,length(sd_banana_fuzzy),1);
            sd_banana(:,1:3) = sd_banana(:,1:3) + repmat(sd_mid_pt+C,length(sd_banana), 1); % placing back in BRIK coord
            %sd_banana_fuzzy(:,1:3) = round(sd_banana_fuzzy(:,1:3));
            sd_banana(:,1:3) = round(sd_banana(:,1:3));

            %[sd_banana] = RBTfm_banana(sd_banana,banana_scaling_factor,tai,p,sd_dist,s_mrkr_pt,sd_banana_mid_pt,s,d,sd_mid_pt,xyz);
              
           % Coloring Banana on the BRIK 
           % This keeps the max intensity under control and increments intensity with banana overlap
           banana_off_screen = 0; % no = 0
           for m=1:length(sd_banana)
               if  (sd_banana(m,1)>mri_dim.x) | (sd_banana(m,2)>mri_dim.y) | (sd_banana(m,3)>mri_dim.z) | (sd_banana(m,1)<1) | (sd_banana(m,2)<1) | (sd_banana(m,3)<1) % Trying to remove any banana edges that go beyond brik dimensions. Dunno why. shouldn't happen. mri_dim.x, mri_dim.y, mri_dim.z.
                   % fprintf('\nError: Source-Deterctor Pair %s-%s element %d  banana point [%d %d %d] exceeds max dimension of [%d %d %d]\n' , char(src(i)), char(det(j)), m, sd_banana(m,1),sd_banana(m,2),sd_banana(m,3),mri_dim.x,mri_dim.y,mri_dim.z);
                   banana_off_screen = 1; %yes=1
               else
                   banana_brik(sd_banana(m,1),sd_banana(m,2),sd_banana(m,3)) =  sd_banana(m,4);
               end
           end
           if banana_off_screen == 1
               fprintf('\nError: Source-Detector Pair %s-%s Off Screen' , src{i}, det{j});
               %fprintf(fid,'\nError: Source-Detector Pair %s-%s Off Screen' , char(src(i)), char(det(j)));
%            elseif banana_off_screen == 0
%                fprintf('\nNo Error: Source-Detector Pair %s-%s' , src{i}, det{j});
%                %fprintf(fid,'\nNo Error: Source-Detector Pair %s-%s' , char(src(i)), char(det(j)));
           end
           
           %for m=1:length(sd_banana_fuzzy)
           %    banana_brik(sd_banana_fuzzy(m,1),sd_banana_fuzzy(m,2),sd_banana_fuzzy(m,3)) = sd_banana_fuzzy(m,4);
           %end
           
           absorbtion_decrement = exp(-mu_a*sd_dist);
           banana_brik = absorbtion_decrement*banana_brik;
           
           % banana_brik = sparse(banana_brik);
           % dc_fact = round(sd_dist/10)
           % banana_brik = decimate(banana_brik,dc_fact);
           % banana_brik = smooth3(banana_brik);
           % banana_brik = interp(banana_brik,dc_fact);
           % banana_brik = full(banana_brik);
           
           % save(['C:\data\mj\' src{i} '-' det{j}], 'banana_brik');
            %pause(0.01);
            
             dig_brik = dig_brik + banana_brik;
             
  %          end %randint if then
             
           
           % Error Checking when Bananas go off screen
%            max_vect = max(sd_banana)
%            if max_vect(1)>mri_dim.x | max_vect(2)>mri_dim.y | max_vect(3)>mri_dim.z
%            %if max(sd_banana(:,1))>mri_dim.x | max(sd_banana(:,2))>mri_dim.y | max(sd_banana(:,3))>mri_dim.z   
%                fprintf('\nError: Source-Deterctor Pair %s-%s' , char(src(i)), char(det(j))) ;
%            else
%                fprintf('\nNo Error: Source-Deterctor Pair %s-%s' , char(src(i)), char(det(j))) ;
%            end
           %scatter3(sd_banana(:,1), sd_banana(:,2), sd_banana(:,3), 2.5, 'r','filled');
           %clear sd_banana;
        else
            %fprintf('\nExcluded Source Detector Pair %s-%s. Distance: %3.0f mm', char(src(i)), char(det(j)), sd_dist);
        end
    end
end

% fprintf('\nScaling Intensities...\n')
if (max(max(max(dig_brik))))==0
    fprintf('banana brik is empty!\n');
else
    dig_brik = dig_brik*max_intensity/(max(max(max(dig_brik))));
end
%fprintf('\nSmoothing...\n');
%dig_brik = smooth3(dig_brik);
if strcmpi(choice.head_mask,'y')
    fprintf('Head Masking...\n');
    dig_brik = dig_brik.*head_mask; % only taking points inside head b/c smoothing gets some outside.
end
% if strcmpi(choice.brain_mask,'y')
%     fprintf('Brain Masking...\n');
%     dig_brik = dig_brik.*brain_mask; % this is to get banana_brik
% end

% dig_brik = permute(dig_brik,[2 1 3]); % done in 'view_bananas'...
dig_brik = floor(dig_brik); % Rounding so can be converted into AFNI brik

% Have something to print and save banana brik as mat file
% choice.image_bananas = input('Image Bananas? (y/n) ', 's');
choice.image_bananas = 'y';
if strcmpi(choice.image_bananas,'y')
    scrsz = get(0,'ScreenSize');
    figure('Position',[scrsz(3)/2 scrsz(4)/2-scrsz(4)/4 ...
                       scrsz(3)/2.2 scrsz(4)/2],'color','w')
    hold on;
    view_bananas(dig_brik,0);
end

banana_time = toc;
% fprintf('\n  Actual Time: %-2.1f min\n', banana_time/60);
% fprintf('\nS-D Bananas Created\n');
% fprintf('\n\nNumber Invalid Channels:  %d\n', (tot_channels - (q-1)));
% fprintf('Number Valid Channels:    %d\n', q-1);
% fprintf('Total Possible Channels:  %d\n', tot_channels);

%  save(['C:\data\mj\label' ], opt_label);
%  save(['C:\data\mj\xyz' ], xyz);
%  save(['C:\data\mj\all_bananas' ], dig_brik);
 
% fprintf('\nBananas Done\n');


